﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace dst
{
    class Matrix<T> where T : Exp
    {
        int row, col;
        T[,] data;

        /// <summary>
        /// 构造一个矩阵
        /// </summary>
        /// <param name="r">行数</param>
        /// <param name="c">列数</param>
        public Matrix(int r, int c)
        {
            row = r;
            col = c;
            data = new T[r, c];
            for (int i = 0; i < row; i++)
            {
                for (int j = 0; j < col; j++)
                    data[i, j] = (T)Exp.Zero;
            }
        }

        /// <summary>
        /// 构造一个方阵
        /// </summary>
        /// <param name="r">阶属</param>
        public Matrix(int r)
            : this(r, r)
        {
        }

        /// <summary>
        /// 方阵的阶数，只读
        /// </summary>
        public int Order
        {
            get
            {
                if (this.row == this.col)
                {
                    return this.row;
                }
                else
                {
                    throw new Exception();
                }
            }
        }

        /// <summary>
        /// 索引器，获取、设置矩阵中的元素
        /// </summary>
        /// <param name="i">所在行</param>
        /// <param name="j">所在列</param>
        /// <returns>元素</returns>
        public T this[int i, int j]
        {
            get { return data[i, j]; }
            set { data[i, j] = value; }
        }

        /// <summary>
        /// 计算矩阵的幂
        /// </summary>
        /// <param name="n">指数</param>
        /// <returns>幂</returns>
        public Matrix<T> Power(int n)
        {
            if (this.row != this.col)
                throw new Exception("Must be square matrix");
            if (n < 0)
                throw new ArgumentException("Exponent must be positive or zero");

            Matrix<T> m = Matrix<T>.Identity(this.Order);

            uint bitmask = 0x80000000;
            while (bitmask > 0)
            {
                if (bitmask > n)
                    continue;

                m = m * m;
                if ((n & bitmask) == 1)
                {
                    m = m * this;
                }
                bitmask >>= 1;
            }

            return m;
        }

        /// <summary>
        /// 计算张量积
        /// </summary>
        /// <param name="b">乘数矩阵</param>
        /// <returns>张量积</returns>
        public Matrix<T> TensorProduct(Matrix<T> b)
        {
            Matrix<T> m = new Matrix<T>(this.row * b.row, this.col * b.col);
            for (int i = 0; i < this.row; i++)
            {
                for (int j = 0; j < this.col; j++)
                {
                    for (int bi = 0; bi < b.row; bi++)
                    {
                        for (int bj = 0; bj < b.col; bj++)
                        {
                            m[i * b.row + bi, j * b.col + bj] = (T)(this[i, j] * b[bi, bj]);
                        }
                    }
                }
            }
            return m;
        }

        /// <summary>
        /// 计算张量积
        /// </summary>
        /// <param name="a">被乘数矩阵</param>
        /// <param name="b">乘数矩阵</param>
        /// <returns>张量积</returns>
        public static Matrix<T> TensorProduct(Matrix<T> a, Matrix<T> b)
        {
            return a.TensorProduct(b);
        }

        /// <summary>
        /// 获得r阶单位矩阵
        /// </summary>
        /// <param name="r">阶数</param>
        /// <returns>单位矩阵</returns>
        public static Matrix<T> Identity(int r)
        {
            Matrix<T> m = new Matrix<T>(r);
            for (int i = 0; i < r; i++)
                m[i, i] = (T)(Exp.Unit);

            return m;
        }

        /// <summary>
        /// 获得r阶上移位矩阵
        /// </summary>
        /// <param name="r">阶数</param>
        /// <returns>上移位矩阵</returns>
        public static Matrix<T> UShift(int r)
        {
            Matrix<T> m = new Matrix<T>(r);
            for (int i = 0; i < r - 1; i++)
                m[i, i + 1] = (T)Exp.Unit;

            return m;
        }

        /// <summary>
        /// 获得r阶n次上移位矩阵
        /// </summary>
        /// <param name="r">阶数</param>
        /// <param name="e">指数</param>
        /// <returns>上移位矩阵</returns>
        public static Matrix<T> UShiftPower(int r, int e)
        {
            Matrix<T> m = new Matrix<T>(r);
            for (int i = 0; i < r - e; i++)
                m[i, i + e] = (T)Exp.Unit;

            return m;
        }

        /// <summary>
        /// 获得r阶下移位矩阵
        /// </summary>
        /// <param name="r">阶数</param>
        /// <returns>下移位矩阵</returns>
        public static Matrix<T> LShift(int r)
        {
            Matrix<T> m = new Matrix<T>(r);
            for (int i = 1; i < r; i++)
                m[i, i - 1] = (T)Exp.Unit;

            return m;
        }

        /// <summary>
        /// 获得r阶e次下移位矩阵
        /// </summary>
        /// <param name="r">阶数</param>
        /// <param name="e">指数</param>
        /// <returns>下移位矩阵</returns>
        public static Matrix<T> LShiftPower(int r, int e)
        {
            Matrix<T> m = new Matrix<T>(r);
            for (int i = e; i < r; i++)
                m[i, i - e] = (T)Exp.Unit;

            return m;
        }

        /// <summary>
        /// 获得r阶副对角线全1的矩阵
        /// </summary>
        /// <param name="r">阶数</param>
        /// <returns>矩阵</returns>
        public static Matrix<T> CIdentity(int r)
        {
            Matrix<T> m = new Matrix<T>(r);
            for (int i = 0; i < r; i++)
                m[i, r - 1 - i] = (T)(Exp.Unit);

            return m;
        }

        public override bool Equals(object obj)
        {
            if ((obj as Matrix<T>) == null)
                return false;
            return Equals((Matrix<T>)obj);
        }

        public bool Equals(Matrix<T> b)
        {
            if (this.row != b.row || this.col != b.col)
                return false;
            for (int i = 0; i < this.row; i++)
                for (int j = 0; j < this.col; j++)
                    if (this[i, j] != b[i, j])
                        return false;
            return true;
        }

        public static bool operator ==(Matrix<T> a, Matrix<T> b)
        {
            if (object.ReferenceEquals(a, null) && object.ReferenceEquals(b, null))
                return true;
            else if (object.ReferenceEquals(a, null) && !object.ReferenceEquals(b, null))
                return false;
            else if (!object.ReferenceEquals(a, null) && object.ReferenceEquals(b, null))
                return false;
            else
                return a.Equals(b);
        }

        public static bool operator !=(Matrix<T> a, Matrix<T> b)
        {
            return !a.Equals(b);
        }

        public static Matrix<T> operator *(Matrix<T> a, Matrix<T> b)
        {
            if (a.col != b.row)
                throw new ArgumentException("can't multiply");

            Matrix<T> m = new Matrix<T>(a.row, b.col);
            for (int i = 0; i < m.row; i++)
            {
                for (int j = 0; j < m.col; j++)
                {
                    m[i, j] = (T)Exp.Zero;
                    for (int k = 0; k < a.col; k++)
                    {
                        if (a[i, k] == Exp.Zero || b[k, j] == Exp.Zero)
                            continue;
                        m[i, j] = (T)(m[i, j] + a[i, k] * b[k, j]);
                    }
                }
            }
            return m;
        }

        public static Matrix<T> operator +(Matrix<T> a, Matrix<T> b)
        {
            if (a.col != b.col || a.row != b.row)
                throw new ArgumentException("can not add");

            Matrix<T> m = new Matrix<T>(a.row, a.col);
            for (int i = 0; i < m.row; i++)
            {
                for (int j = 0; j < m.col; j++)
                {
                    m[i, j] = (T)(a[i, j] + b[i, j]);
                }
            }
            return m;
        }

        public override string ToString()
        {
            string s = "";
            StringBuilder sb = new StringBuilder();
            int[] width = new int[this.col];
            string[,] tmpM = new string[this.row, this.col];

            for (int i = 0; i < this.row; i++)
            {
                for (int j = 0; j < this.col; j++)
                {
                    string tmp = this[i, j].ToString();
                    tmpM[i, j] = tmp;
                    if (width[j] < tmp.Length)
                        width[j] = tmp.Length;
                }
            }
            for (int i = 0; i < this.row; i++)
            {
                for (int j = 0; j < this.col; j++)
                {
                    sb.Append(string.Format("{0,-" + width[j] + "}", tmpM[i, j]));
                    if (j != (this.col - 1))
                    {
                        sb.Append(", ");
                    }
                }
                sb.Append("\n");
            }
            s = sb.ToString();
            return s;
        }

        /// <summary>
        /// 判断本矩阵是否为Toeplitz矩阵
        /// </summary>
        /// <returns></returns>
        public bool isToeplitz()
        {
            for (int i = 1; i < this.row; i++)
            {
                for (int j = 1; j < this.col; j++)
                {
                    if (this[i, j] != this[i - 1, j - 1])
                    {
#if DEBUG
                        Console.WriteLine(String.Format("Matrix[{0}, {1}] != Matrix[{2}, {3}]", i, j, i - 1, j - 1));
#endif
                        return false;
                    }
                }
            }
            return true;
        }

        /// <summary>
        /// 如果确定乘积是Toeplitz矩阵，可以用这个函数加速乘法
        /// </summary>
        /// <param name="b"></param>
        /// <returns></returns>
        public Matrix<T> Mul2Toep(Matrix<T> b)
        {
            Matrix<T> r = new Matrix<T>(this.row);

            for (int i = 0; i < r.col; i++)
            {
                for (int k = 0; k < r.row; k++)
                {
                    if (this[0, k] == Exp.Zero || b[k, i] == Exp.Zero)
                        continue;
                    r[0, i] = (T)(r[0, i] + this[0, k] * b[k, i]);
                }

                int j = i + 1;
                while (j < r.col)
                {
                    r[j - i, j] = r[0, i];
                    j++;
                }
            }

            for (int i = 1; i < r.row; i++)
            {
                for (int k = 0; k < r.row; k++)
                {
                    if (this[i, k] == Exp.Zero || b[k, 0] == Exp.Zero)
                        continue;
                    r[i, 0] = (T)(r[i,0] + this[i, k] * b[k, 0]);
                }

                int j = i + 1;
                while (j < r.col)
                {
                    r[j, j - i] = r[i, 0];
                    j++;
                }
            }
            return r;
        }

        public override int GetHashCode()
        {
            return base.GetHashCode();
        }
    }
}
